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AN  EFFICIENT  NUMERICAL  SOLUTION  TECHNIQUE  FOR  WAVE  PROPAGATION 
IN  HORIZONTALLY  STRATIFIEO  OCEAN  ENVIRONMENTS 

hy 

Henrik  Schmidt  and  Finn  B.  Jensen 


ABSTRACT 


A  new  solution  technique  for  wave  propagation  in  horizontally  stratified 
viscoelastic  media  is  presented.  The  model  provides  a  full -wave  solution 
for  the  field  generated  by  a  single  source,  as  well  as  for  that  generated 
by  a  vertical  source  array.  It  allows  the  spatial  distribution  of  the 
acoustic  field  to  be  evaluated  at  least  one  order  of  magnitude  faster  than 
with  existing  models  based  on  the  Thomson-Haskel 1  solution  technique. 
Computational  examples  demonstrate  the  model's  versatility  in  providing 
exact  solutions  to  a  wide  range  of  guided  and  non-guided  propagation 
problems  in  underwater  acoustics. 


INTRODUCTION 


A  numerical  model  providing  a  full -wave  solution  to  propagation  problems  in 
horizontally  stratified  fluid/solid  environments  is  presented.  During  the 
last  decade  a  number  of  similar  models  have  been  developed,  known  in  under¬ 
water  acoustics  as  Fast  Field  Programs  <l-3>  and  in  seismology  as  full  wave 
field  response  models  <4>.  These  models  are,  however,  in  general  based  on 
the  Thomson-Haskel 1  matrix  method,  which  allows  for  field  calculations  for 
only  one  source/recei ver  combination  at  a  time.  Hence  calculations  of  the 
field  as  a  function  of  depth  and  of  fields  produced  by  vertical  source 

arrays  need  several  separate  runs,  with  the  calculation  time  being  propor¬ 
tional  to  the  number  of  sources  and  receivers. 

In  this  model  a  more  direct  approach  is  taken.  The  field  is  considered  as 
a  superposition  of  two  fields:  one  produced  by  the  sources  in  the  absence 
of  boundaries,  and  an  unknown  field  satisfying  the  homogeneous  wave 

equation.  The  unknown  field  Is  then  determined  from  the  boundary  con¬ 

ditions  to  be  satisfied  at  each  Interface.  This  approach  Is  equivalent  to 
the  technique  used  by  Ewing  et  al  <5>  to  derive  analytic  solutions, 

although  not  In  closed  form,  for  special  cases  with  only  a  few  layers.  The 
model  described  here  is  In  fact  just  a  general  numerical  Implementation  of 
their  approach,  but  compared  with  the  Thomson-Haskel 1  technique  there  are  a 
number  of  Important  advantages.  First  there  are  no  restrictions  on  the 
number  of  sources,  as  the  fields  produced  by  more  sources  within  a  layer 
are  simply  superimposed.  Second,  any  number  of  receiver  depths  can  be 
treated  with  one  solution,  since  the  unknown  field  is  found  In  all  layers 
simultaneously.  Thus  the  spatial  distribution  of  the  field  produced  by  a 
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single  source  or  by  a  beam-generating  vertical  array  can  be  determined  In  a 
computationally  efficient  manner.  Third,  even  for  only  one  source/receiver 
combination,  this  solution  technique  as  Implemented  here  yields  a  code  that 
Is  an  order  of  magnitude  faster  than  existing  models  using  the 
Thomson -Haskell  technique.  This  is  partly  due  to  the  fact  that  the 
Thomson-Haskel  1  technique  requires  additional  specialized  numerical  effort 
in  order  to  ensure  stability  in  cases  where  evanescent  waves  are  propa¬ 
gating  in  very  thick  layers.  Such  problems,  by  contrast,  are  removed  auto¬ 
matically  from  this  direct  global  solution  technique  by  choosing  a  proper 
local  coordinate  system  within  each  layer. 

We  describe  the  theoretical  basis  for  wave  propagation  in  stratified 
viscoelastic  media  in  Sect.  1,  following  closely  Ewing  et  al  <5>.  The 
numerical -solution  technique  particular  to  this  model  is  described  In 
Sect.  2,  with  details  on  the  numerical  stability  being  addressed  In 
Appendix  A.  Section  3  is  dedicated  to  computational  examples  for  three 
important  propagation  problems  in  underwater  acoustics:  1)  Lloyd-mlrror 
Interference  effects  at  the  sea  surface,  2)  low-frequency  propagation  In  a 
shallow-water  duct,  and  3)  reflection  of  a  narrow  beam  of  sound  at  the  sea 
floor.  The  paper  ends  with  a  summary  and  conclusions. 


DERIVATION  OF  FIELD  EQUATIONS 


The  representation  of  the  field  in  terms  of  integral  solutions  closely 
follows  the  presentation  given  by  Ewing  et  al  <5>  and  thus  only  an  outline 
is  given  here.  The  environment  is  assumed  to  be  horizontally  stratified 
and  all  layers,  including  the  upper  and  lower  half  spaces,  are  considered 
to  be  isotropic  and  homogeneous  viscoelastic  continua  with  Lam§  constants 
Xn  and  pn  and  density  pn.  The  subscript  n  refers  to  the  layer  number. 

The  field  equations  are  here  derived  in  cylindrical  geometry.  The  deriva¬ 
tion  of  the  corresponding  field  representations  in  plane  geometry  is  given 
elsewhere  <6>. 

A  cylindrical  coordinate  system  {r,e,z}  is  introduced,  with  the  z-axis  per¬ 
pendicular  to  the  interfaces  and  positive  downwards,  see  Fig.  1. 
In  the  absence  of  body  forces  the  equations  of  motion  will  be  satisfied 
if  the  displacement  components  {u,v,w}  in  layer  'n'  are  expressed  in  terms 
of  the  scalar  potentials  {®n,  ¥n,  An}  as 
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FIG.  1  HORIZONTALLY  STRATIFIED  MEDIUM. 
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If  the  medium  Is  a  fluid,  pn  vanishes,  and  only  the  potential  *n  Is  pre¬ 
sent,  This  can  however  be  considered  as  a  special  case  of  the  general 
solid  case,  and  only  when  necessary  Is  It  treated  separately  In  the 
following. 

The  sources  are  assumed  to  be  harmonic  and  vibrate  with  the  common  angular 
frequency  u>  .  In  complex  notation  a  time  dependence  of  the  form  exp(lut) 
Is  assumed,  a  common  factor  that  will  be  omitted  in  the  following.  The 
viscoelastic  attenuation  can  be  accounted  for  by  letting  the  Lam£  constants 
be  complex.  The  wave  equations,  Eqs.  4  and  5,  now  take  the  form 


(y2+h^  ®n  =  0  (Eq.  ®) 

{*n.  An}  -  0  ,  (Eq.  9) 


where  hn  and  kn  are  the  wavenumbers  for  compressional  and  shear  waves  respec¬ 
tively: 


2 
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A„+2pn 


kn  *  (C~)2 
n  Csn 
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nn 


(Eq.  10) 

(Eq.  11) 


If  only  compressional  sources  placed  on  the  axis  r=0  are  allowed,  the  field 
will  be  axlsymmetric  with  vanishing  tangential  displacements  v(r,z).  The 
torsional  potential  ¥n  can  then  be  excluded  from  the  analysis. 

By  applying  the  Hankel  transform  to  Eqs.  8  and  9  the  following  Integral 
representations  are  obtained  for  the  solutions  <5>: 


*n(r  >*) 

•7  [A’(s)e'Zan(s) 

+  A+(s)eZan^S^]  sJ0(rs)ds 

(Eq.  12) 
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(Eq.  14) 
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Note  that  Eqs.  12  and  13  are  just  decompositions  of  the  total  fields  In  up 
and  downgoing  conical  waves  integrated  over  all  angles  (wavenumbers). 

Because  of  the  attenuation  and  the  choice  of  time  factor,  the  Imaginary 
parts  of  hn  and  kn  are  negative,  and  the  complex  s-plane  Is  cut  along 
branches  going  from  hn  and  kn  to  -i«  and  from  -hn  and  -kn  to  +i«. 

Substituting  Eqs.  12  and  13  into  Eqs.  1  and  3  gives  the  following  Integral 
representations  for  the  displacements: 


u(r.z)  |  _  *  /  (-s  A'  e'Zan  -  s  A+  eZa" 
o  n  n 

+  $n  B”  e  zen  -  en  B*  eZ6n}  sJ^rsJds  (Eq.  16) 

w(r,2)  )  n  *  l  (-«n  e'Zan  +  an  A+  eZan 

+  sB“  e‘ZBn  +  s  B+  eZPn}  sJ0(rs)ds.  (Eq.  17) 


The  stress  components  involved  in  the  boundary  conditions  follow  by 
Hooke's  law: 
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(Eq.  18) 
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-(2s2-k2)(B;  e’Z0n  +  B+  eZ6n)}  sJj (rs)ds. (Eq.  19) 


In  the  case  of  a  fluid  layer  the  displacements  follow  directly  from  Eqs,  16 

and  17  by  setting  B"  and  B+  to  zero.  The  shear  stress  opt  vanishes, 

n  n 

whereas  Eq.  18  has  to  be  replaced  by 
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(Eq.  20) 


To  obtain  expressions  for  the  total  field  in  layer  V,  the  field 
produced  by  the  sources  within  the  layer  has  to  be  added.  Only 
congressional  sources  at  r=0  are  considered  and  the  field  produced  by 
each,  in  an  infinite  medium  with  the  material  properties  of  layer  V,  has 
I  the  integral  representation  <5>: 


-I  z’zs  “n 


*  00  p  ' 

*n(r,z)  =  /  - 

o  an 


sJ0(rs)ds, 


An(r,z)  =  0  , 


(Eq.  21) 

(Eq.  22) 


where  zs  is  the  source  depth. 

The  corresponding  displacements  involved  in  the  boundary  conditions  are 
again  obtained  from  Eqs.  1  and  3  as: 


u(r,z) 
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an 


sJi (rs)ds 


(Eq.  23) 
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(Eq.  24) 


while  the  stresses  are  given  by  Hooke's  law: 

■  |  z_zs  |  °n 

®zz(r*z)  I  n  =  un  I  (2s2 -k2)  - -  sJ0(rs)ds  (Eq.  25) 

bn  an 

°xz(p.z)  |  n  =  un  /  2s  s1gn(z-zs)e"  ^  Z'Zs  ^  °n  sJ^rsJds.  (Eq.  26) 
o 

In  the  fluid  case,  Eq.  25  is  replaced  by 
*  •  “  I  z-zs  I  “n 

ozz(r.z)  |  n  *  Xn  hn  /  ® — — -  sJ0(rs)ds.  (Eq.  27) 


If  more  than  one  source  Is  present  In  the  layer,  the  kernels  In  Eqs.  23  to 
27  are  replaced  by  a  sum  over  the  number  of  sources. 
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The  field  at  each  Interface  now  has  two  different  integral  representations, 
one  from  the  layer  above  and  one  from  the  layer  below.  The  boundary  con¬ 
ditions  have  to  be  satisfied  at  all  ranges  r;  thus  they  have  to  be 
satisfied  by  the  kernels  in  the  integral  representations  as  well,  leading 
to  a  linear  system  of  equations  in  the  unknown  arbitrary  functions  A+,  A", 
B+  and  B".  This  system  can  of  course  be  solved  analytically,  as  done  for 
special  cases  In  <5>,  but  for  a  general  multilayered  environment  a  numeri¬ 
cal  solution  is  more  convenient.  This  numerical  solution  and  the  sub¬ 
sequent  evaluation  of  the  integral  transforms  is  treated  in  the  next 
section. 


2  NUMERICAL  SOLUTION  TECHNIQUE 

The  numerical  solution  is  divided  into  two  parts.  First  the  unknown 
arbitrary  functions  are  found  at  a  discrete  number  of  horizontal  wavenum¬ 
bers  from  the  system  of  equations  that  expresses  the  boundary  conditions  to 
be  satisfied.  Secondly  the  field  is  found  at  selected  depths  by  evaluation 
of  the  Integral  transforms.  The  first  part  is  the  most  critical  in  rela¬ 
tion  to  computation  time,  and  this  is  the  part  where  the  present  model  dif¬ 
fers  In  approach  from  earlier  models  of  the  same  type.  Most  emphasis  is 
therefore  given  to  this  part  in  the  following. 

The  layers  are  numbered  from  1  to  N  (Fig.  1)  where  the  upper  halfspace  is 
number  1  and  the  lower  halfspace  is  number  N.  If  the  kernels  for  the  field 
parameters  Involved  in  the  boundary  conditions  at  interface  'm',  and 
evaluated  in  layer  V,  are  expressed  in  vector  form  as  (v}m,  then  the 

local  system  of  equations  to  be  satisfied  at  interface  'm'  is 


{v}m  -  {vr  ♦  ivr  *  {*}■♦  m* 

m  m  m+i  m+l 


(Eq.  28) 


where  the  asterisk  denotes  the  source  contribution.  Here  and  in  the 
following  the  subscript  refers  to  the  layer  number  and  the  superscript  to 
the  Interface  number.  In  the  general  solid/solid  case,  Eq.  28  expresses 
continuity  of  u,  w,  o2z  and  orz.  If  one  of  the  media  is  a  fluid,  w,  ozz  and 
o rz  has  to  be  continuous,  whereas  the  fluid/fluid  case  requires  only  con¬ 
tinuity  of  w  and  czz.  Where  one  of  the  media  is  a  vacuum  the  stresses  ozz 
and  a rz  must  vanish.  The  number  of  elements  in  the  vectors  (v}m  there¬ 
fore  vary  from  Interface  to  Interface.  n 

The  contributions  from  the  unknown  homogeneous  solutions  are  now  isolated 
on  the  left-hand  side,  and  Eq.  28  can  be  rewritten  as 


(Eq.  29) 


where  [c]n  is  a  coefficient  matrix  and  {a}  is  a  vector  including  the 
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degrees  of  freedom  for  layer  'm',  i.e.  all  or  some  of  the  arbitrary 

functions  A-,  A+,  B"  and  B+. 

m  m  m  m 

To  combine  the  N-l  local  systems  of  equations  (one  for  each  Interface)  Into 
one  global  system,  the  following  unique  mapping  is  introduced: 

{a}m  =  [S]m  {A}  (Eq.  30) 

{V}  «  V  [T]1  {v}1  .  (Eq.  31) 

i=l 


Here  {A}  and  {V}  are  global  vectors  containing  all  degrees  of  freedom  and 
field  parameters  involved  in  the  boundary  conditions.  The  matrices  [S]m 
and  [Tji  are  extremely  sparse,  containing  only  zeroes  and  ones.  These 
matrices  are  very  similar  to  the  topology  matrices  known  from  the  finite- 
element  technique,  and  as  is  the  case  there,  they  are  never  set  up  in  the 
actual  computer  code,  but  replaced  by  a  set  of  pointers  <7>. 

Insertion  of  the  mappings  of  Eqs.  30  and  31  into  Eq.  29  yields  the 
following  global  system  of  equations 


[ C] { A }  »  (V}  ,  (Eq.  32) 


where  the  coefficient  matrix  is 


N-l 


(Cl  -  l  t[T]'[c]][S).- 
i  =1 


(Eq.  33) 


The  right-hand  side  due  to  the  source  contributions  is 

*  N-l  *  .  *  . 

{V}  -  l  [T]i  ({v}!  -  {v}1)  . 

i =1  i+l  i 


(Eq.  34) 


As  the  mappings  of  Eqs.  30  and  31  are  unique,  the  summations  and  the  matrix 
multiplications  in  Eqs.  33  and  34  never  need  be  performed,  but  can  be 
replaced  by  a  unique  set  of  pointers  connecting  the  elements  of  the  global 
system  with  those  of  the  similar  local  systems,  as  illustrated  in  Fig.  2. 

The  pointers  defined  by  the  mappings  of  Eqs.  30,  31,  33  and  34  are  depen¬ 
dent  only  on  the  environment  and  can  thus  be  determined  a  priori .  The 
calculations  needed  at  each  horizontal  wavenumber  are  then  delimited  to  the 
creation  of  the  local  coefficient  matrices  and  right-hand  sides,  and  to  the 
solution  of  Eq.  32.  Then  the  kernels  In  Eqs.  16  to  20  and  23  to  27  can  be 
determined  for  any  number  of  depths. 
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The  local  coefficient  matrices  are  functions  of  the  horizontal  wavenumber 
s  and  the  depth  z  of  the  actual  interface.  As  can  be  observed  from 
Eqs.  16  to  20,  the  dependence  on  z  Is  present  only  in  the  exponential 
function;  thus  the  two  local  matrices  for  layer  number  a  can  be  written 
as 


[c]*-1  =  [d]  [e]*-l 

t  it 


(Eq.  35) 


[el*  =  [d]  [ej* 
£  £  £ 


(Eq.  36) 


where  [d]*,  is.  a  matrix  that  is  a  function  of  the  horizontal  wavenumber 

only,  and  [e]i  is  a  diagonal  matrix  containing  the  exponential  functions 
£ 

for  interface  number  i  evaluated  in  layer  number  £.  By  introducing  a 
local  coordinate  system  where  the  plane  z=0  coincides  with  one  of  the 
interfaces,  the  number  of  necessary  calculations  of  the  complex  exponential 
function  can  be  reduced.  If,  for  example,  the  interface  above  the  layer  is 
chosen,  Eqs.  35  and  36  take  the  simpler  forms: 

[c]*-1  =  [d],  (Eq.  37) 


EcJf  *  [d]t[e]‘  .  (Eq.  38) 


To  build  the  global  coefficient  matrix  by  means  of  the  pointers  defined  by 
Eq.  33,  it  is  necessary  only  to  calculate  for  each  layer  the  elements  of 
the  matrix  [d]* .which  are  very  simple  functions  of  the  horizontal  wave- 
number,  and  the  exponentials  in  [e]|  . 

Although,  from  a  theoretical  point  of  view,  the  origin  of  the  local  coor¬ 
dinate  system  can  be  chosen  arbitrarily,  its  choice  is  quite  important  for 
the  numerical  stability  when  large  real  arguments  of  the  exponential  func¬ 
tions  appear.  This  is  so  for  thick  layers  and  large  values  of  the  horizon¬ 
tal  wavenumber.  In  these  cases  a  numerical  solution  based  on  the  original 
Thomson-Haskel  1  technique  becomes  unstable,  as  is  clear  from  the  physical 
significance  of  the  exponential  functions.  These  express  the  depth  beha¬ 
viour  of  the  field,  and  when  the  arguments  are  Imaginary  they  correspond  to 
up  and  downgoing  conical  waves.  However,  when  the  arguments  become  real, 
as  is  the  case  for  large  horizontal  wavenumbers,  the  waves  are  Inhomoge¬ 
neous  conical  waves  travelling  In  the  horizontal  direction  with  an  exponen¬ 
tial  decay  away  from  the  Interfaces  (provided  that  no  sources  are  present 
In  the  layer).  When  the  arguments  become  very  large,  the  field  produced  at 
the  opposite  Interface  by  these  exponential  "tails"  will  vanish.  At  these 
wavenumbers  a  thick  layer  will  therefore  behave  exactly  as  an  Infinite 
halfspace.  As  the  Thomson-Haskel 1  technique  requires  that  the  solution  Is 
"propagated"  through  all  layers,  numerical  stability  problems  obviously 
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arise  because  this  is  impossible  at  these  wavenumbers.  By  factorizing  out 
the  exponential  functions  a  stable  solution  can  be  obtained,  and  several 
computational  schemes  for  this  have  been  proposed,  see  for  instance 
Franssens  <8>.  These  modifications  of  the  original  technique,  however, 
yield  a  more  complex  and  time-consuming  computer  code. 

With  the  present  technique  it  is  possible  to  obtain  unconditionally  stable 
solutions  by  applying  gaussian  elimination  with  partial  pivoting  to  Eq.  32 
and  by  choosing  proper  origins  of  the  local  coordinate  systems  together 
with  a  proper  mapping.  This  is  demonstrated  in  detail  in  Appendix  A. 

The  fact  that  no  special  numerical  effort,  except  the  pivoting,  is  needed 
to  ensure  numerical  stability  in  this  technique  is  one  of  the  major  reasons 
for  the  order-of -magnitude  improvement  in  computational  speed  over  that  of 
the  modified  Thomson -Haskel 1  technique  <2>.  Another  important  factor  is 

the  mapping  technique,  which  ensures  that  only  a  few  simple  functions  of 

the  horizontal  wavenumber  have  to  be  calculated  for  each  layer.  Due  to  the 
special  band  structure  of  the  global  coefficient  matrix  (Fig.  2),  gaussian 
elimination  with  partial  pivoting  can  be  performed  very  efficiently. 

An  important  advantage  of  the  present  technique  is  that  most  operations  can 

be  easily  vectorized,  making  it  very  well  suited  for  implementation  on  an 

array  processor.  Because  the  local  coefficient  matrices,  Eqs.  37  and  38, 
are  similar  for  all  layers,  the  calculation  of  the  elements,  including  the 
square  roots  and  the  exponentials,  can  be  vectorized.  The  indexed  move 
operations  of  the  mapping  are  carried  out  very  efficiently  on  array  pro¬ 
cessors,  and  the  gaussian  elimination  is  by  nature  a  sequence  of  vector 
operations.  The  implementation  of  this  technique  on  an  FPS-164  array  pro¬ 
cessor  has  therefore  yielded  a  considerable  improvement  in  computational 
speed  compared  with  serial  computers,  as  can  be  observed  from  the  calcula¬ 
tion  times  given  in  Sect.  3. 

When  the  system  of  equations  (Eq.  32)  has  been  solved,  the  kernels  In 
Eqs.  16  to  20  and  23  to  27  can  be  evaluated  at  any  depth,  z,  with  the  only 
additional  functions  needed  being  exp(±z  an);  all  other  functions  and 
expressions  have  been  evaluated  while  setting  up  the  system  of  equations. 
The  present  technique  is  therefore  very  efficient  if  the  field  Is  to  be 
determined  at  many  different  depths. 

It  is  well  known  that  for  guided  propagation  in  loss-less  media  the  kernels 
In  the  integral  representations  will  have  poles  on  the  real  axis 
corresponding  to  normal  modes,  Rayleigh  waves,  etc.  In  these  cases  a 
direct  numerical  Integration  along  the  real  axis  is  inconvenient.  There 
are  two  ways  In  which  this  problem  can  be  overcome.  One  is  to  deform  the 
contour  of  integration  into  the  complex  plane;  the  other  is  to  Introduce 
physically  realistic  attenuations  in  the  layers,  which  will  result  In  the 
poles  moving  out  into  the  complex  plane  and  consequently  make  real-axis 
Integration  possible.  The  latter  approach  has  the  advantage  that  the 
Integration  can  be  performed  by  means  of  the  Fast  Field  Technique  <1>,  and 
this  is  the  choice  made  here.  It  should  be  stressed  however  that  the  solu¬ 
tion  technique  described  here  is  valid  for  both  complex  or  real  values  of 
the  horizontal  wavenumber  's',  and  that  any  complex  integration  contour 
could  be  chosen. 


ll 
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The  Bessel  functions  are  expressed  in  terms  of  Hankel  functions  as 

J*(rs)  =I[H(l)(rs)  +  H(2)(rs)]  (Eq.  39) 

and  each  integral  is  split  into  .two.  As  only  outgoing  waves  are  con¬ 
sidered,  the  integrals  involving  H(1)(rs)  are  neglected,  and  H(z)(rs)  Is 
replaced  by  its  asymptotic  form  *  * 


H(2)(rs)  ~ 

*  rs+» 


-ICrs-d*!)!] 
e  2  2 


(Eq.  40) 


As  shown  in  <3>  this  approximation  yields  insignificant  errors  for  ranges 
beyond  a  few  wavelengths. 

The  field  integrals  now  take  the  form 

F(r,z)  =  -L-  ei(^7  J  f(s,z).  /s  e‘irsds  ,  (Eq.  41) 

a/STF  o 

where  4=0  for  w  and  o22  and  £=1  for  u  and  or2.  The  integral  in  Eq.  41  can 
be  evaluated  by  means  of  an  FFT,  but  in  order  to  do  this,  the  interval  of 
integration  must  be  finite.  However,  because  of  the  exponential  decay  of 
the  kernels  for  s  going  towards  infinity,  the  truncation  error  can  be  made 
arbitrarily  small.  The  truncated  wavenumber  space  is  discretized  as 
follows 


sn  =  s0  +  nAs  ,  n  =  0,  1,...  (M-l). 


(Eq.  42) 


In  addition,  the  range  interval  of  interest  is  discretized  as: 


rm  =  r0  =  m&r  ,  m  =  0,  1,  ...  (M-l) 


(Eq.  43) 


where 

Aras  =  —  (Eq.  44) 

M 

and  M  is  a  power  of  2.  The  field  integrals  are  now  evaluated  as 
F(rm,2)  AS  g-i [Sn^m" (t+l/2)ff /2] 

*Y  [f (sn»z)e"1r°nAs  /s^]  e'i2nmn/M  ,  (Eq.  45) 

n=0 


where  the  summation  is  performed  by  means  of  an  FFT  yielding  the  field  at 
all  M  ranges,  Eq.  43,  simultaneously. 
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In  the  case  of  plane  geometry  the  approximation  given  by  Eq.  40  Is  not 
needed,  and  the  FFT  integration  can  be  performed  directly  <6>. 


3  COMPUTATIONAL  EXAMPLES 

In  this  section  we  present  numerical  solutions  to  three  Important  propaga¬ 
tion  problems  In  underwater  acoustics:  1)  the  Lloyd-mlrror  Interference 
effect  at  the  sea  surface,  2)  low-frequency  propagation  In  a  shallow-water 
duct,  and  3)  reflection  of  a  narrow  beam  of  sound  at  the  sea  floor.  These 
problems  have  been  selected  so  as  to  demonstrate  the  full  capability  of  the 
numerical  model  described  earlier,  i.e.  its  applicability  to  a  wide  range 
of  guided  and  non-guided  propagation  problems  In  both  plane  and  cylindrical 
geometry.  The  computational  efficiency  of  the  code  will  become  apparent 
from  the  quoted  calculation  times  on  both  a  VAX-11/750  and  on  a  FPS-164 
array  processor. 


3.1  Sea-surface  Interference  effects 


This  is  a  classical  problem  dealt  with  In  most  text  books  on  underwater 
acoustics  <9>.  We  consider  an  omnidirectional  line  source  In  a  halfspace 
limited  above  by  a  perfectly  reflecting  boundary.  The  acoustic  field  In 
the  halfspace  (Fig.  3a)  becomes  quite  complex,  exhibiting  a  series  of 
Interference  beams  radiating  out  at  different  angles  (In  the  contour  plot 
black  Indicates  high  sound  intensity  and  white  low  sound  Intensity).  A 
closed-form  solution  to  this  problem  can  be  easily  derived  <9>,  showing 
that  the  beam  pattern  arises  as  an  interference  effect  between  the  two 
possible  sound  paths  from  source  to  receiver,  namely  the  direct  path  and 
the  surface-reflected  path.  This  interference  pattern  Is  often  referred 
to  as  Lloyd-mirror  beams,  with  the  beam  directions  far  from  the  source 
given  by 


sin  0m  =  (2m-l)  A/4zs 


where  A  Is  the  acoustic  wavelength  and  zs  the  source  depth;  6m  Is  the 
beam  angle  relative  to  horizontal. 

The  numerical  solution  given  in  Fig.  3a  is  in  plane  geometry  with  a  line 
source  placed  33.75  m  below  the  sea  surface.  The  frequency  Is  100  Hz  and 
the  water  sound  speed  is  1500  m/s.  The  source  depth  Is  2.25  A,  leading  to 
five  Lloyd-mlrror  beams  at  angles  given  by  sin  9m  =  (2m-l)/9  or  6m  *  6.4, 
19.5,  33.7,  51.1,  and  90.0°.  The  computed  angular  spectrum  for  this  propa¬ 
gation  problem  Is  shown  in  Fig.  3b,  and  we  see  that  energy  is  In  fact  pro¬ 
pagating  In  the  directions  given  by  the  simple  formula,  and  that  the  peak 
Intensity  Is  the  same  in  all  five  beams.  The  beamwldth,  however.  Increases 
with  Increasing  beam  angle. 

The  Iso-Intensity  contour  plot  of  Fig.  3a  is  an  exact  numerical  solution 
to  the  sea-surface  reflection  problem.  For  a  spatial  grid  of  50  x  1024 
points  the  computation  time  was  5  min  on  the  VAX-11/750  and  18  s  on  the 
FPS-164. 
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Even  though  the  field  solution  here  is  available  in  closed  form,  many  of 
the  standard  wave-theory  models  in  ocean  acoustics  are  unable  to  handle 
this  particular  problem.  This  is  so  for  all  normal -mode  models  that  do  not 
include  the  branch-line  integral  (continuous  spectrum)  contribution  to  the 
field  solution.  It  Is  also  so  for  all  parabolic  equation  models,  with 
their  Inherent  small-angle  approximation  limiting  the  solution  validity  to 
at  best  ±  40°  away  from  the  horizontal.  This  would  here  mean  that  only  the 
first  three  Lloyd-mirror  beams  would  be  correctly  resolved  with  the  parabo¬ 
lic  equation.  In  fact  we  consider  this  particular  test  problem  extremely 
well-suited  for  checking  parabolic  equation  solutions,  since  inherent 
amplitude  and  phase  errors  affect  both  beam  directions  and  beam  inten¬ 
sities. 


3.2  Propagation  in  a  shallow-water  duct 

An  isovelocity  shallow-water  duct  that  is  limited  above  by  the  sea  surface 
and  below  by  a  homogeneous  penetrable  bottom  constitutes  a  waveguide  with 
Interesting  frequency-dependent  propagation  characteristics.  We  will  solve 
this  problem  in  cylindrical  geometry,  with  environmental  parameters  as 
follows:  100  m  water  depth  with  a  speed  of  1500  m/s,  a  solid  bottom  with 
compressional  speed  of  1800  m/s,  shear  speed  of  600  m/s,  compresslonal 
attenuation  of  0.1  dB/x,  and  shear  attenuation  of  0.2  dB/x.  The  density 
ratio  between  bottom  and  water  is  2.0.  We  have  intentional ly  used  a  solid 
bottom  in  order  to  demonstrate  that  the  often-neglected  shear  properties 
can  be  of  considerable  importance  in  low-frequency  ocean  acoustics. 

A  ful  1 -spectrum  solution  at  30  Hz  is  displayed  in  Fig.  4  for  a  source  at 
95  m  and  a  receiver  at  100  m  (on  the  bottom).  The  propagating  energy  ver¬ 
sus  horizontal  wavenumber  is  shown  on  the  integrand  plot  of  Fig.  4a,  which 
displays  the  kernel  of  Eq.  12  with  source  contributions  included.  We  have 
here  Identified  three  spectral  regions,  defined  in  terms  of  the  phase  velo¬ 
cities  as: 


Evanescent  spectrum:  0  <  cph  <  1500  m/s 

Discrete  spectrum  :  1500  <  cp^  <  1800  m/s 
Continuous  spectrum:  1800  <  cph  , 

where  the  phase  velocity  Cnh  is  related  to  the  horizontal  wavenumber  's' 
through  Cph  -  u/s,  with  w  Deing  the  angular  frequency  of  the  source.  In 
terms  of  an  angular  spectrum  the  evanescent  waves  (modes)  have  non-real 
angles  at  the  receiver  depth.  The  discrete  modes,  on  the  other  hand, 
correspond  to  propagation  angles  between  horizontal  (cph=1500  m/s)  and 
33.6°(cph=1800  m/s),  while  the  continuous  modes  correspond  to  steep  propa¬ 
gation  angles  above  33.6°.  Thus  it  is  the  compressional  wave  speeds  In 
water  and  bottom  that  determine  the  spectral  regions. 

We  see  from  the  Integrand  (Fig.  4a)  that  there  are  five  peaks  in  the 
spectrum  corresponding  to  five  preferred  modes  of  propagation.  There  are 
two  discrete  modes,  two  continuous  (virtual)  modes,  and  one  evanescent 
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mode,  which  can  be  shown  to  be  an  Interface  wave  of  the  Scholte  type 

<10,11,12>  propagating  along  the  sea  floor.  The  interface  wave  has  a  phase 
velocity  of  530  m/s,  which  is  slower  than  any  of  the  body  waves  of  this 

problem  (1500,  1800,  600  m/s).  Moreover  it  has  its  highest  amplitude  at 
the  water/bottom  interface,  with  exponentially  decaying  amplitude  away  from 
the  interface.  The  existence  of  the  interface  wave  is  intrinsically 
related  to  the  shear  properties  of  the  bottom,  and,  hence,  if  we  neglect 
shear  we  automatically  exclude  this  propagation  path. 

The  integrand  plot  (Fig.  4a)  shows  that  the  two  discrete  modes  have  highest 
excitation  and  lowest  loss  (poles  located  close  to  the  real  axis  in  the 
complex  wavenumber  space).  The  discrete  modes  will  therefore  dominate  pro¬ 
pagation  beyond  the  nearfield,  as  shown  in  the  lower  transmission-loss 

curve  (shear  »  600  m/s)  of  Fig.  4b.  Beyond  2  km  we  have  a  typical  2-mode 
interference  pattern  caused  by  the  discrete  modes.  The  two  highly  atte¬ 
nuated  virtual  modes  are  responsible  for  the  irregularities  of  the  inter¬ 
ference  pattern  at  short  ranges,  while  the  interface  mode  causes  the 

“high-frequency"  noise  on  the  propagation-loss  curve.  The  above  detailed 
analysis  of  the  contributions  of  the  various  spectral  components  to  the 
propagation-loss  curve  can  be  demonstrated  explicitly  by  solving  for  just 
selected  parts  of  the  wavenumber  spectrum. 

A  full-spectrum  calculation  without  shear  shows  better  propagation  at  long 
ranges,  but  with  essentially  the  same  interference  structure  (Fig.  4b). 

Hence  the  effect  of  shear  is  here  to  increase  losses  for  waterborne  propa¬ 
gation. 

We  now  move  to  a  frequency  of  5  Hz  (Fig.  5),  which  is  below  the  cut-off 
frequency  for  discrete  modes.  Here  the  effect  of  shear  is  to  make  propaga¬ 
tion  conditions  better  than  they  are  in  the  case  where  shear  is  neglected 
(see  Fig.  5b).  Get  us  shortly  analyze  the  spectral  content  of  the 

integrand  plot  (Fig.  5a).  There  are  no  discrete  modes  and  only  one  highly 
damped  virtual  mode.  The  important  peak,  however,  is  associated  with  the 
interface  wave,  which  has  low  attenuation  and  therefore  becomes  the  most 
Important  propagation  path  at  this  frequency.  The  resulting  propagation 
loss  (Fig.  5b),  with  a  shear  speed  of  600  m/s,  shows  interference  between 
the  virtual  mode  and  the  interface  mode.  However,  since  the  virtual  mode 
is  highly  attenuated,  only  the  interface  mode  remains  at  long  ranges 
(beyond  5  km).  If  we  neglect  shear  there  is  no  interface  mode;  propagation 
thus  becomes  very  poor  because  energy  can  then  propagate  only  in  the  con¬ 
tinuous  mode.  On  the  other  hand,  we  notice  from  Figs.  4  and  5  that  when 

shear  is  included,  propagation  is  just  as  good  at  5  Hz  as  it  is  at  30  Hz, 
even  though  the  propagation  mechanisms  are  quite  different  in  the  two 
cases.  The  importance  of  shear  in  low-frequency  ocean  acoustics  has  been 
addressed  in  more  detail  elsewhere  <10,11>. 

The  numerical  model  described  in  this  paper  is  extremely  well  suited  for 
solving  low-frequency  guided  wave-propagation  problems.  A  ful 1 -spectrum 
solution  is  readily  available  with  minimum  computational  effort.  A 
propagation-loss  curve  ,  such  as  shown  in  Fig.  4  or  5,  is  produced  in  65  s 
on  the  VAX-11/750  and  in  just  5  s  on  the  FPS-164.  These  computation  times 
relate  to  the  use  of  2048  integration  points  along  the  wavenumber  axis. 

The  most  commonly  used  models  for  solving  guided  propagation  problems  are 
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FIG.  5  HORIZONTAL  WAVENUMBER  SPECTRUM  FOR  SHALLOW  WATER  CASE  AT 
5  Hz  (a),  AND  ASSOCIATED  TRANSMISSION  LOSS  CURVES  (b). 
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based  on  normal -mode  theory.  Again,  if  such  a  model  does  not  include  the 
branch-line  integral,  the  continuous  spectrum  is  not  included.  Moreover, 
most  normal-mode  models  solve  a  real  eigenvalue  problem,  and  hence  neglect 
shear  altogether,  or  at  best  include  shear  in  a  perturbational  manner  as 
just  a  loss  mechanism.  In  any  event,  the  interface  wave  is  not  part  of  the 
modal  solution,  and  hence  such  a  model  would  fail  when  applied  to  a  propa¬ 
gation  problem  not  dominated  by  discrete  modes. 

The  ful 1 -spectrum  solution  technique  presented  here  is  clearly  superior  to 
most  existing  models  in  handling  guided-wave  propagation  in  horizontally 
stratified  fluid/solid  layers.  Only  in  the  case  of  many  guided  low-loss 
modes  do  numerical  problems  arise,  since  many  integration  points  are  then 
needed  in  order  to  resolve  the  peaks  in  the  integrand. 


3.3  Beam  reflection  at  a  water/bottom  interface 

In  this  last  example  we  study  the  reflection  and  transmission  of  a  narrow 
beam  of  sound  at  a  water/bottom  interface.  This  example  has  been  chosen  to 
demonstrate  that  directional  sources  (beams)  are  handled  by  the  numerical 
model  in  a  computationally  efficient  manner.  The  beam  is  generated  by  a 
vertical  source  array  composed  of  a  number  of  equidistantly  spaced  line 
sources  (plane  geometry),  and  the  resulting  acoustic  field  is  found  by 
superposition  of  the  contributions  from  the  individual  sources,  as 
described  in  Sect.  1.  The  beam  direction  is  varied  by  appropriately 
phasing  the  source  elements,  and  the  intensity  distribution  across  the  beam 
can  be  selected  by  applying  an  amplitude  weighting  across  the  array.  By 
varying  the  array  distance  from  the  interface  and  the  number  of  source  ele¬ 
ments  (half-wavelength  spacing),  a  beam  of  arbitrary  width  can  be 
generated.  Moreover,  beams  can  be  focused  or  defocused  by  phasing  the 
sources  appropriately.  By  using  a  "physical"  array  of  finite  length  to 
generate  the  beam  we  automatical ly  ensure  that  only  realistic  beams  are 
generated,  i.e.  a  focused  beam  in  water  (with  realistic  attenuation)  cannot 
be  infinitely  narrow,  because  that  would  lead  to  infinite  intensity  in  the 
focal  plane,  which  is  unphysical.  Hence  this  model  treats  an  entirely 
realistic  physical  system. 

The  interest  in  the  beam-reflection  problem  goes  back  to  an  experimental 
study  by  Muir  et  al  <13>,  where  it  was  shown  in  a  laboratory  experiment 
that  a  narrow  beam  of  sound  impinging  on  the  bottom  at  grazing  angles  below 
the  critical  angle  will  not  be  totally  reflected  as  predicted  for  an  infi¬ 
nitely  wide  plane  wave.  Instead  there  is  significant  energy  transmitted 
into  the  bottom,  this  energy  increasing  as  the  incident  beam  narrows.  This 
interesting  observation  has  remained  unexplained  for  years,  but  we  shall 
demonstrate  numerically  that  the  beam  penetration  is  a  simple  consequence 
of  the  angular  spectrum  shape  associated  with  narrow  beams. 

The  beam  reflection/transmission  problem  will  be  solved  in  plane  geometry 
for  an  environment  similar  to  the  one  used  by  Muir  et  al  <13>.  We  consider 
a  water  halfspace  overlying  a  fluid-bottom  halfspace.  The  sound  speeds  are 
1450  m/s  in  the  water  and  1675  m/s  in  the  bottom,  giving  a  critical 
(grazing)  angle  of  exactly  30°.  The  density  ratio  between  bottom  and  water 
is  2.0  and  we  initially  neglect  attenuation  in  the  bottom.  Calculations 
are  done  for  a  frequency  of  20  kHz,  which  gives  an  acoustic  wavelength  (x) 
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In  the  water  of  7.25  cm.  We  have  chosen  a  source  array  of  121  elements 
with  A/2  spacing,  giving  an  array  length  of  4.25  m.  A  gausslan  amplitude 
weighting  was  applied  across  the  array,  producing  a  gausslan  beam  with  low 
si delobe  levels. 

Calculated  field  solutions  are  given  In  Fig.  6  for  a  focused  beam  that  Is 
2.5  A  wide  at  the  Interface.  The  beamwldth  Is  measured  across  the  beam 
between  the  3-dB  down  points.  The  arrows  In  the  upper  part  of  Fig.  6  Indi¬ 
cate  the  beam  directions;  also  shown  are  the  angles  of  Incidence  (6j), 
reflection  (er),  and  transmission  (6*),  all  measured  along  the  direction  of 
maximum  amplitude  In  a  beam  with  respect  to  horizontal.  Three  cases  are 
considered,  e-f  =  25.0°,  30.0°,  35.0°,  with  the  critical  angle  being  30°. 
Hence  the  upper  beam  Is  Incident  at  5°  below  the  critical  angle,  and  a 
transmitted  beam  Is  still  present  In  the  bottom  (0 t  *  11.5°).  With 
increasing  angle  of  incidence  more  energy  Is  propagated  into  the  bottom 
(contours  show  losses  In  arbitrary  dB's).  We  also  notice  that  the  angle  of 
reflection  is  lower  than  the  angle  of  incidence  (non-specular  reflection). 
These  beam  reflection  and  transmission  properties  are  In  qualitative 
agreement  with  the  experimental  results  reported  by  Muir  et  al  <13>,  even 
though  the  experiment  was  carried  out  with  a  true  3-dlmenslonal  pencil 
beam,  while  the  calculations  here  are  done  for  a  2-dlmenslonal  beam. 

A  full  understanding  of  the  observed  beam  behaviour  can  be  obtained  by  com¬ 
paring  the  spectral  content  of  the  Incident  beam  with  the  angular  reflec¬ 
tivity  characteristics  of  the  Interface,  as  shown  in  Fig.  7.  The  upper 
figure  displays  angular  spectra  for  beams  of  different  widths  and  Incident 
at  the  critical  angle.  Note  that  the  narrow  beam  has  the  widest  spectrum, 
while  the  spectral  width  decreases  with  Increasing  beamwldth.  In  this 
angular  representation  an  infinitely  wide  plane  wave  becomes  a  delta  func¬ 
tion  with  only  one  direction  of  propagation.  Any  beam  of  finite  width  has 
a  spectrum  of  finite  width;  the  very  narrow  beams  have  broad  spectra.  Indi¬ 
cating  that  energy  Is  propagating  over  a  wide  range  of  angles.  Hence  It  Is 
a  contradiction  to  state  that  narrow  beams  are  highly  directional  <13,14>. 

Let  us  concentrate  on  the  wide  spectrum  In  Fig.  7a  (D*  2.5  a),  which  has 
Its  peak  energy  propagating  at  30°  but  also  has  significant  energy  propa¬ 
gating  at  20°  and  40°.  Hence  many  propagation  directions  are  associated 
with  the  narrow  beam.  Figure  7b  shows  the  plane-wave  reflection  loss  and 
phase  shift  versus  angle  at  the  Interface,  indicating  perfect  reflection 
below  30°  with  a  phase  shift,  and  an  Increasing  reflection  loss  above  30° 
with  no  phase  shift.  Comparing  Figs.  7a  and  7b  shows  that  the  left 
(hatched)  half  of  the  beam  spectra  are  perfectly  reflected  (though  with  a 
phase  shift)  while  the  right  half  Is  partly  transmitted.  Thus  In  practice 
we  need  to  weight  the  Incident -beam  spectrum  with  the  reflectivity  charac¬ 
teristics  of  the  Interface  In  order  to  get  the  spectrum  for  the  reflected 
or  transmitted  beam.  The  above  simple  argument  explains  qualitatively  what 
to  expect.  For  example,  the  hatched  part  of  the  spectrum  will  be  entirely 
reflected,  while  the  right  half  will  be  partly  transmitted.  The  different 
weightings  put  on  the  two  half-spectra  mean  that  the  peak  energy  In  the 
reflected  beam  moves  to  smaller  angles  (0r  <  0 1 ) .  as  also  observed  in 
Fig.  4.  On  the  other  hand,  only  energy  In  the  right  half  of  the  spectrum 
will  be  transmitted,  and  the  beam  direction  (6t)  can  be  determined  from  the 
peak  In  the  transmitted  spectrum. 
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FIG.  6  REFLECTION  AND  TRANSMISSION  OF  NARROW  BEAM  OF  SOUND  (D~2.5\)  AT 

WATER/ BOTTOM  INTERFACE.  THE  CRITICAL  ANGLE  IS  3CP  AND  THERE  IS  NO 
ATTENUATION  IN  THE  BOTTOM.  'S'  S'' 
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It  is  clear  from  Fig.  7  that  even  though  each  spectral  component  is 
reflected  and  transmitted  according  to  Snell's  law,  the  total  beam,  in 

which  phase  and  amplitude  weightings  are  applied  to  the  broad  spectrum, 

cannot  be  expected  to  be  so.  On  the  other  hand,  the  beam  results  should 
approach  Snell's  law  for  increasing  beamwidth  (decreasing  spectral  width). 
In  fact,  as  shown  in  Fig.  8,  we  do  retrieve  Snell's  law  in  the  limit  of 

very  wide  beams. 

Finally,  when  a  realistic  bottom  loss  of  0.8  dB/x  is  included  we  notice  no 
substantial  effect  on  the  calculated  beam  angles,  even  though  the 
transmitted  beam  becomes  strongly  attenuated. 

One  last  point  will  be  briefly  addressed.  It  was  noted  in  the  experimental 
results  <13>  that  the  transmitted  beam  was  horizontally  displaced,  i.e.  the 
incident  and  transmitted  beam  axes  intersect  the  interface  at  different 
points.  A  closer  look  at  Fig.  6  shows  that  both  the  reflected  and  the 
transmitted  beam  centres  are  displaced;  this  is  because  a  lateral  wave  is 
excited  when  energy  is  incident  on  the  interface  below  the  critical  angle. 
The  reflected  and  transmitted  fields  are  then  composed  of  contributions 
from  both  the  "specular"  beam  and  the  lateral  wave  field,  causing  an 
apparent  lateral  displacement  of  the  beams.  The  displacement  can  be  both 
forward  and  backward  and  will  be  a  function  of  beamwidth  as  well  as  of 
angle  of  incidence.  Furthermore,  reflected  and  transmitted  beams  will 
generally  have  different  displacements. 

The  computational  efficiency  of  the  numerical  code  is  particularly  evident 
for  this  case,  in  which  field  calculations  from  121  sources  have  been  made 
on  a  spatial  grid  of  50  x  512  points.  Each  contoured  beam  plot  in  Fig.  6 
takes  6  min  on  the  VAX-11/750  and  21  s  on  the  FPS-164.  Not  only  are  these 
acceptable  calculation  times,  but  this  code  seems  to  be  the  only  one 
capable  of  doing  these  types  of  calculations. 

The  results  presented  here  are  exact  numerical  solutions  for  2-dimensional 
beams,  against  which  approximate  solutions  <14>  can  be  checked.  Similar 
beam  reflection  studies  at  fluid/solid  interfaces  have  been  reported  by  the 
authors  elsewhere  <6,15,16>. 


SUMMARY  AND  CONCLUSIONS 

We  have  presented  an  efficient  numerical  solution  technique  for  general 
applications  to  wave  propagation  in  horizontally  stratified  viscoelastic 
media.  The  model  provides  a  full -wave  solution  for  the  field  generated  by 
a  single  source  as  well  as  for  that  generated  by  a  vertical  source  array. 
It  allows  the  spatial  distribution  of  the  acoustic  field  to  be  evaluated  at 
least  one  order  of  magnitude  faster  than  with  existing  models  based  on  the 
Thomson-Haskel  1  solution  technique.  The  computational  examples  clearly 
demonstrate  the  model's  versatility  In  providing  exact  solutions  to  a  wide 
range  of  guided  and  non-guided  propagation  problems  in  underwater 
acoustics,  but  the  model  is  equally  applicable  to  problems  in  ultrasonics 
and  selsmlcs.  Moreover,  the  computational  efficiency  of  the  numerical  code 
has  made  It  feasible  (on  array  processors)  to  perform  pulse  calculations 
based  on  Fourier  synthesis  of  time  signals,  even  when  hundreds  of  frequency 
samples  are  needed. 


Grazing  angle  (deg) 


FIG.  8  DEVIATION  FROM  SNELL' S  LAM  FOR  REFLECTION  AND  TRANSMISSION  OF 

NARROW  BEAM  OF  SOUND  AT  WATER/BOTTOM  INTERFACE.  Bottom  attenuation 
is  neglected  and  the  indicated  beamwidths  are  measured  between  the 
3  dB  down  points  of  the  incident  beam  at  the  interface. 
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APPENDIX  A 
NUMERICAL  STABILITY 


The  physically  realistic  attenuation  Introduced  in  the  environment  makes 
the  global  system  of  equations  (Eq.  32)  theoretically  well  conditioned  for 
real  horizontal  wavenumbers.  Because  of  the  limited  accuracy  of  digital 
computers,  this  does  not,  however,  guarantee  stability  of  the  numerical 
solution  technique.  In  this  paper  gaussian  elimination  is  used  to  solve 
Eq.  32,  which,  in  general,  requires  scaling  and  pivoting  to  ensure  numeri¬ 
cal  stability. 

Equations  16  to  20  of  the  main  text  show  that  the  difference  in  dimension 
between  the  displacements  and  the  stresses  can  yield  several  orders  of 
magnitude  In  difference  between  the  coefficients  in  the  corresponding 
Eqs.  29  and  32.  In  such  cases  gaussian  elimination  with  partial  pivoting 
will  not  ensure  unconditional  stability,  owing  to  truncation  errors.  The 
equations  therefore  have  to  be  scaled  In  order  to  make  all  coefficients 
dimensionless.  The  choice  of  actual  scaling  factors  is  not  critical.  Here 
all  stress  equations  are  divided  by  the  constant  u>  ps,  where  ps  is  the 
density  of  any  iayer.  The  displacement  equations  are  made  dimensionless  by 
dividing  by  the  wavenumber  hs  for  the  same  layer. 

For  purely  imaginary  arguments  of  the  exponential  functions,  the  magnitude 
of  the  exponentials  is  unity,  and  the  solution  will  be  unconditionally 
stable.  For  real  arguments,  i.e.  the  evanescent  regime  with  s>hn  or  kn, 
large  differences  in  order  of  magnitude  can  however  appear  again.  In  this 
case  the  differences  appear  within  the  rows,  and  scaling  of  rows  is  there¬ 
fore  not  applicable.  Scaling  of  columns  could  theoretically  be  applied, 
but  this  requires  considerable  additional  computations  at  each  horizontal 
wavenumber,  and  the  problem  of  exponential  overflow  in  the  actual  computer 
would  still  remain.  With  a  little  help  from  the  physical  significance  of 
the  problem  It  is  however  straightforward  to  obtain  stability  also  In  this 
case,  without  any  additional  computational  effort. 

This  Is  most  easily  demonstrated  by  means  of  an  example.  Consider  a  series 
of  fluid  layers  with  a  very  thick  layer  V  below  the  lowermost  source.  In 
the  evanescent  case,  I.e.  s»hn,  this  layer  will  behave  like  an  infinite 

half  space,  which  means  that  the  arbitrary  function  A+  associated  with  the 

n 

exponentially  Increasing  depth  function  (eZan)  should  vanish.  The  solution 
of  the  global  system  of  equations  {Eq.  32)  by  means  of  gaussian  elimination 
with  partial  pivoting  can  be  forced  to  yield  this  behaviour  automatically 
simply  by  choosing  a  mapping  that  yields  the  structure  of  the  system  showed 
In  Fig.  Al,  The  upper  boundary  of  the  layer,  I.e.  Interface  n-1.  Is 
selected  as  the  origin  of  the  local  coordinate  system,  which  removes  the 
exponentials  from  the  submatrix  [c]nn~l.  When  the  argument  to  the  exponen¬ 
tial  functions  In  [c]n  becomes  large  It  is  simply  truncated;  the  column 
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Fig.  A.l  MAPPING  ENSURING  NUMERICAL  STABILITY  IN  EVANESCENT  CASE 
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closest  to  the  diagonal  of  the  global  system  then  contains  at  least  one 
very  large  number,  whereas  the  lefthand  column  contains  two  almost 
vanishing  numbers.  No  pivoting  will  then  be  performed  between  the  sub- 

matrices  [c]n_*  and  [c]n  .  As  no  sources  are  present  in  the  layers  >n,  the 
n  n 

right-hand  sides  corresponding  to  interfaces  >n  will  then  always  be  zero. 

The  element  of  [c]n  placed  on  the  diagonal  of  the  global  system  remains 
n 

very  large  (eventually  after  pivoting  with  the  row  below)  compared  with  the 
other  elements  In  that  row  and  will  therefore  force  the  corresponding 

arbitrary  function  A+  to  be  zero  after  triangulation  and  back-substitution. 

n 

The  lowest  equation  for  Interface  n-1  then  yields  the  correct  value  of 

A-,  and  back-substitution  can  continue, 
n 

In  a  similar  way  it  can  be  shown  that  for  layers  above  the  uppermost  source 
the  origin  of  the  local  coordinate  system  must  be  placed  at  the  lower  boun¬ 
dary  in  order  to  ensure  numerical  stability  if  the  columns  of  the  sub¬ 
matrices  are  organized  as  shown  in  Fig.  Al.  Stability  problems  can  then 
appear  only  in  two  cases.  The  first  is  when  sources  are  placed  within  a 
very  thick  layer,  but  this  can  be  avoided  by  adding  dummy  interfaces  above 
and  below  the  sources.  The  other  is  the  rather  unrealistic  situation  In 
which  sources  are  present  on  both  sides  of  a  thick  layer.  The  solution  of 
such  a  problem  would  require  two  separated  solutions  if  the  evanescent 
regime  Is  considered. 

In  the  case  of  solid  layers,  unconditional  stability  is  ensured  in  the  same 
way  but  because  two  large  exponentials  can  appear  in  any  layer  It  Is 
slightly  more  complicated  to  demonstrate. 
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